******************************************************************************************************************************************
*** REPLICATION FILE FOR APPENDIX "ON THE ECONOMIC CONSEQUENCES OF MASS SHOOTINGS" BY ABEL BRODEUR AND HASIN YOUSAF
*** PREPARED FOR PUBLICATION IN THE REVIEW OF ECONOMICS AND STATISTICS
******************************************************************************************************************************************

clear all
set more off
set mat 3000
global path = "Change path here"

******************************************************************************************************************************************
*** TABLE A1 
******************************************************************************************************************************************

*** COLUMN 1
use "$path\finaldata.dta", clear

sort fips year

foreach var of varlist ln_ca1_pop_1 ln_emp_pop ln_real_qp1_pop ln_est_pop ln_real_qp1_job ln_job_est lnhpi pop25g white black hispanic lhs coll {
reg L1.ch`var' mse if aroundsuccfailms==1 & postattack==0
}

*** COLUMN 2
foreach var of varlist ln_ca1_pop_1 ln_emp_pop ln_real_qp1_pop ln_est_pop ln_real_qp1_job ln_job_est lnhpi pop25g white black hispanic lhs coll {
reg L1.ch`var' mse if aroundmsneighbor==1 & postattack==0
}

*** COLUMN 3
sort fips year
foreach var of varlist ln_emp_pop ln_real_qp1_pop ln_est_pop ln_ca1_pop_1 lnVIOLENT lnPROPERTY male hs racediv poverty gini lnhom_gun lnsui_gun {
g lag`var'=L1.`var'
}

global lag_dependentvar = "lagln_emp_pop lagln_real_qp1_pop lagln_est_pop"
global lag_controls_crime = "laglnVIOLENT laglnPROPERTY"
global lag_controls_census = "lagmale laghs lagracediv lagpoverty laggini"
global lag_controls_gun = "laglnhom_gun laglnsui_gun"
global controls_geo = "capital_state coastal_county large_hub medium_hub"

***
probit ms lagln_ca1_pop_1 $lag_dependentvar $lag_controls_crime $lag_controls_gun i.year, cluster(fips)
predict p_score, pr
sort p_score year lagln_emp_pop lagln_real_qp1_pop lagln_est_pop fips

g matchms = 0
replace matchms = 1 if ms[_n+1]==1 & mse!=1 
replace matchms = 1 if ms[_n+2]==1 & mse!=1 
replace matchms = 1 if ms[_n+3]==1 & mse!=1 

sort fips year
by fips: egen matchmse = max(matchms)

g aroundmatchms = 0
forvalues j=-6(1)6{
by fips: replace aroundmatchms=1 if ms[_n+`j']==1
by fips: replace aroundmatchms=1 if matchms[_n+`j']==1
}

reg L1.chln_ca1_pop_1 mse if aroundmatchms==1 & postattack==0
outreg2 using "balancingtest5", replace excel
foreach var of varlist ln_emp_pop ln_real_qp1_pop ln_est_pop ln_real_qp1_job ln_job_est lnhpi pop25g white black hispanic lhs coll {
reg L1.ch`var' mse if aroundmatchms==1 & postattack==0
outreg2 using "balancingtest5", append excel
}

******************************************************************************************************************************************
*** TABLE A2
******************************************************************************************************************************************

use "$path\finaldata.dta", clear

reghdfe ln_emp_pop successful if aroundms==1, absorb(fips i.div_9_all##i.year) cluster(fips)
g esample=(e(sample)==1)
keep if esample==1

reghdfe ln_real_qp1_job successful if aroundms==1, absorb(fips i.div_9_all##i.year) cluster(fips)
reghdfe ln_job_est successful if aroundms==1, absorb(fips i.div_9_all##i.year) cluster(fips)

******************************************************************************************************************************************
*** TABLE A3
******************************************************************************************************************************************

use "$path\finaldata.dta", clear

reghdfe ln_emp_pop successful if aroundms==1, absorb(fips i.div_9_all##i.year) cluster(fips)
g esample=(e(sample)==1)
keep if esample==1

reghdfe ln_emp_pop_qcew_p successful if aroundms==1, absorb(fips i.div_9_all##i.year) cluster(fips)
reghdfe ln_emp_pop_qcew_gf successful if aroundms==1, absorb(fips i.div_9_all##i.year) cluster(fips)
reghdfe ln_emp_pop_qcew_edu successful if aroundms==1, absorb(fips i.div_9_all##i.year) cluster(fips)
reghdfe ln_emp_pop_qcew_oth successful if aroundms==1, absorb(fips i.div_9_all##i.year) cluster(fips)

******************************************************************************************************************************************
*** TABLE A4
******************************************************************************************************************************************

use "$path\finaldata.dta", clear
sort fips year

reghdfe ln_emp_pop successful13-successful17 successful19-successful25 if aroundms==1, absorb(fips i.div_9_all##i.year) cluster(fips)

reghdfe ln_real_qp1_pop successful13-successful17 successful19-successful25 if aroundms==1, absorb(fips i.div_9_all##i.year) cluster(fips)

reghdfe ln_est_pop successful13-successful17 successful19-successful25 if aroundms==1, absorb(fips i.div_9_all##i.year) cluster(fips)

******************************************************************************************************************************************
*** TABLE A5
******************************************************************************************************************************************
use "$path\finaldata.dta", clear

keep if aroundms==1

did_multiplegt ln_emp_pop fips year successful, breps(1000) placebo(5) dynamic(5) robust_dynamic 
did_multiplegt ln_real_qp1_pop fips year successful, breps(1000) placebo(5) dynamic(5) robust_dynamic 
did_multiplegt ln_est_pop fips year successful, breps(1000) placebo(5) dynamic(5) robust_dynamic

******************************************************************************************************************************************
*** TABLE A6
******************************************************************************************************************************************

use "$path\finaldata.dta", clear

egen countyid=group(fips)
egen yeart=group(year) 

g e=1
bysort fips: egen e1=sum(e)
drop if e1!=21 & e1!=24 // To ensure it is a balanced panel
drop e e1
xtset countyid yeart, delta(1)
keep fips year countyid yeart ln_emp_pop ln_real_qp1_pop ln_est_pop successful
drop if ln_emp_pop==.
compress

bacondecomp ln_emp_pop successful, ddetail
bacondecomp ln_real_qp1_pop successful, ddetail
bacondecomp ln_est_pop successful, ddetail

******************************************************************************************************************************************
*** TABLE A7
******************************************************************************************************************************************
use "$path\finaldata.dta", clear

reghdfe ln_emp_pop successful if aroundms==1, absorb(fips year i.state_all##i.year) cluster(fips)
reghdfe ln_real_qp1_pop successful if aroundms==1, absorb(fips i.state_all##i.year) cluster(fips)
reghdfe ln_est_pop successful if aroundms==1, absorb(fips i.state_all##i.year) cluster(fips)

******************************************************************************************************************************************
*** TABLE A8
******************************************************************************************************************************************

*** - PANEL A
use "$path\finaldata.dta", clear

reghdfe ln_emp_pop successful postattack if aroundsuccfailms==1, absorb(fips i.state_all##i.year) cluster(fips)
reghdfe ln_real_qp1_pop successful postattack if aroundsuccfailms==1, absorb(fips i.state_all##i.year) cluster(fips)
reghdfe ln_est_pop successful postattack if aroundsuccfailms==1, absorb(fips i.state_all##i.year) cluster(fips)

*** - PANEL B
use "$path\finaldata.dta", clear

reghdfe ln_emp_pop successful if aroundmsneighbor==1, absorb(fips i.state_all##i.year) cluster(fips)
reghdfe ln_real_qp1_pop successful if aroundmsneighbor==1, absorb(fips i.state_all##i.year) cluster(fips)
reghdfe ln_est_pop successful if aroundmsneighbor==1, absorb(fips i.state_all##i.year) cluster(fips)

*** - PANEL C
use "$path\finaldata.dta", clear

sort fips year
foreach var of varlist ln_emp_pop ln_real_qp1_pop ln_est_pop ln_ca1_pop_1 lnVIOLENT lnPROPERTY male hs racediv poverty gini lnhom_gun lnsui_gun {
g lag`var'=L1.`var'
}

global controls = "lagln_emp_pop lagln_real_qp1_pop lagln_est_pop laglnVIOLENT laglnPROPERTY lagmale laghs lagracediv lagpoverty laggini laglnhom_gun laglnsui_gun capital_state coastal_county large_hub medium_hub"

***
probit ms lagln_ca1_pop_1 $controls i.year, cluster(fips)
predict p_score, pr
sort p_score year lagln_emp_pop lagln_real_qp1_pop lagln_est_pop fips

g matchms = 0
replace matchms = 1 if ms[_n+1]==1 & mse!=1 
replace matchms = 1 if ms[_n+2]==1 & mse!=1 
replace matchms = 1 if ms[_n+3]==1 & mse!=1 

sort fips year
by fips: egen matchmse = max(matchms)

g aroundmatchms = 0
forvalues j=-6(1)6{
by fips: replace aroundmatchms=1 if ms[_n+`j']==1
by fips: replace aroundmatchms=1 if matchms[_n+`j']==1
}

reghdfe ln_emp_pop successful if aroundmatchms==1, absorb(fips i.state_all##i.year) cluster(fips)
reghdfe ln_real_qp1_pop successful if aroundmatchms==1, absorb(fips i.state_all##i.year) cluster(fips)
reghdfe ln_est_pop successful if aroundmatchms==1, absorb(fips i.state_all##i.year) cluster(fips)

******************************************************************************************************************************************
*** TABLE A9
******************************************************************************************************************************************
use "$path\finaldata.dta", clear

reghdfe ln_emp_pop successful if aroundms==1 & msn<=1, absorb(fips year) cluster(fips)
reghdfe ln_emp_pop successful if aroundms==1 & msn<=1, absorb(fips i.region_4_all##i.year) cluster(fips)
reghdfe ln_emp_pop successful if aroundms==1 & msn<=1, absorb(fips i.div_9_all##i.year) cluster(fips)

reghdfe ln_real_qp1_pop successful if aroundms==1 & msn<=1, absorb(fips year) cluster(fips)
reghdfe ln_real_qp1_pop successful if aroundms==1 & msn<=1, absorb(fips i.region_4_all##i.year) cluster(fips)
reghdfe ln_real_qp1_pop successful if aroundms==1 & msn<=1, absorb(fips i.div_9_all##i.year) cluster(fips)

reghdfe ln_est_pop successful if aroundms==1 & msn<=1, absorb(fips year) cluster(fips)
reghdfe ln_est_pop successful if aroundms==1 & msn<=1, absorb(fips i.region_4_all##i.year) cluster(fips)
reghdfe ln_est_pop successful if aroundms==1 & msn<=1, absorb(fips i.div_9_all##i.year) cluster(fips)

******************************************************************************************************************************************
*** TABLE A10
******************************************************************************************************************************************

use "$path\finaldata.dta", clear

*** PANEL A
reghdfe ln_emp_pop successful if aroundms==1 & msustoday==0, absorb(fips year) cluster(fips)
reghdfe ln_emp_pop successful if aroundms==1 & msustoday==0, absorb(fips i.region_4_all##i.year) cluster(fips)
reghdfe ln_emp_pop successful if aroundms==1 & msustoday==0, absorb(fips i.div_9_all##i.year) cluster(fips)

reghdfe ln_real_qp1_pop successful if aroundms==1 & msustoday==0, absorb(fips year) cluster(fips)
reghdfe ln_real_qp1_pop successful if aroundms==1 & msustoday==0, absorb(fips i.region_4_all##i.year) cluster(fips)
reghdfe ln_real_qp1_pop successful if aroundms==1 & msustoday==0, absorb(fips i.div_9_all##i.year) cluster(fips)

reghdfe ln_est_pop successful if aroundms==1 & msustoday==0, absorb(fips year) cluster(fips)
reghdfe ln_est_pop successful if aroundms==1 & msustoday==0, absorb(fips i.region_4_all##i.year) cluster(fips)
reghdfe ln_est_pop successful if aroundms==1 & msustoday==0, absorb(fips i.div_9_all##i.year) cluster(fips)

*** PANEL B
g aroundms2 = 0
replace aroundms2=1 if msustoday==1
replace aroundms2=1 if mse==0

reghdfe ln_emp_pop successful if aroundms2==1, absorb(fips year) cluster(fips)
reghdfe ln_emp_pop successful if aroundms2==1, absorb(fips i.region_4_all##i.year) cluster(fips)
reghdfe ln_emp_pop successful if aroundms2==1, absorb(fips i.div_9_all##i.year) cluster(fips)

reghdfe ln_real_qp1_pop successful if aroundms2==1, absorb(fips year) cluster(fips)
reghdfe ln_real_qp1_pop successful if aroundms2==1, absorb(fips i.region_4_all##i.year) cluster(fips)
reghdfe ln_real_qp1_pop successful if aroundms2==1, absorb(fips i.div_9_all##i.year) cluster(fips)

reghdfe ln_est_pop successful if aroundms2==1, absorb(fips year) cluster(fips)
reghdfe ln_est_pop successful if aroundms2==1, absorb(fips i.region_4_all##i.year) cluster(fips)
reghdfe ln_est_pop successful if aroundms2==1, absorb(fips i.div_9_all##i.year) cluster(fips)


******************************************************************************************************************************************
*** TABLE A11
******************************************************************************************************************************************

*** PANEL A
use "$path\finaldata.dta", clear
merge 1:1 fips year using "$path\ms2.dta" 

sort fips year
by fips: egen ms2n = sum(ms2)
by fips: egen ms2e = max(ms2)

by fips: egen minyr2 = min(year) if ms2==1
by fips: egen minyr = min(minyr2) if ms2e==1
by fips: g yrms2 = year if ms2==1
by fips: replace yrms2 = yrms2[_n-1] if yrms==.
by fips: replace yrms2 = minyr if yrms2==. & ms2n!=0 & year<=yrms2
g yearsince2 = year - yrms2
g post2 = (yearsince2>=0 & yearsince2!=.)
drop minyr2 minyr
g mspost2 = ms2e*post2

sort fips year
g aroundms2 = 0
forvalues j=-6(1)6{
by fips: replace aroundms2=1 if ms2[_n+`j']==1
}
replace aroundms2=1 if ms2e==0

reghdfe ln_emp_pop mspost2 if aroundms2==1, absorb(fips i.div_9_all##i.year) cluster(fips)
reghdfe ln_real_qp1_pop mspost2 if aroundms2==1, absorb(fips i.div_9_all##i.year) cluster(fips state_all)
reghdfe ln_est_pop mspost2 if aroundms2==1, absorb(fips i.div_9_all##i.year) cluster(fips)


*** PANEL B
use "$path\finaldata.dta", clear
merge 1:1 fips year using "$path\ms3.dta" 

sort fips year
by fips: egen ms3n = sum(ms3)
by fips: egen ms3e = max(ms3)

by fips: egen minyr2 = min(year) if ms3==1
by fips: egen minyr = min(minyr2) if ms3e==1
by fips: g yrms3 = year if ms3==1
by fips: replace yrms3 = yrms3[_n-1] if yrms3==.
by fips: replace yrms3 = minyr if yrms3==. & ms3n!=0 & year<=yrms3
g yearsince3 = year - yrms3
g post3 = (yearsince3>=0 & yearsince3!=.)
drop minyr2 minyr
g mspost3 = ms3e*post3

sort fips year
g aroundms3 = 0
forvalues j=-6(1)6{
by fips: replace aroundms3=1 if ms3[_n+`j']==1
}
replace aroundms3=1 if ms3e==0

reghdfe ln_emp_pop mspost3 if aroundms3==1, absorb(fips i.div_9_all##i.year) cluster(fips)
reghdfe ln_real_qp1_pop mspost3 if aroundms3==1, absorb(fips i.div_9_all##i.year) cluster(fips)
reghdfe ln_est_pop mspost3 if aroundms3==1, absorb(fips i.div_9_all##i.year) cluster(fips)

*** PANEL C
use "$path\finaldata.dta", clear

drop if succdeaths==4
sort fips year
g aroundms2 = 0
forvalues j=-6(1)6{
by fips: replace aroundms2=1 if ms[_n+`j']==1
}
replace aroundms2=1 if mse==0

reghdfe ln_emp_pop mspost if aroundms2==1, absorb(fips i.div_9_all##i.year) cluster(fips)
reghdfe ln_real_qp1_pop mspost if aroundms2==1, absorb(fips i.div_9_all##i.year) cluster(fips state_all)
reghdfe ln_est_pop mspost if aroundms2==1, absorb(fips i.div_9_all##i.year) cluster(fips)

*** PANEL D
use "$path\finaldata.dta", clear

drop if succdeaths==4 | succdeaths==5
sort fips year
g aroundms2 = 0
forvalues j=-6(1)6{
by fips: replace aroundms2=1 if ms[_n+`j']==1
}
replace aroundms2=1 if mse==0

reghdfe ln_emp_pop mspost if aroundms2==1, absorb(fips i.div_9_all##i.year) cluster(fips)
reghdfe ln_real_qp1_pop mspost if aroundms2==1, absorb(fips i.div_9_all##i.year) cluster(fips state_all)
reghdfe ln_est_pop mspost if aroundms2==1, absorb(fips i.div_9_all##i.year) cluster(fips)

************************************************************************************************************************
*** TABLE A12
************************************************************************************************************************

use "$path\finaldata.dta", clear

reghdfe ln_emp_pop successful postattack sage if aroundsuccfailms==1, absorb(fips i.div_9_all##i.year smale sweapon) cluster(fips)
reghdfe ln_real_qp1_pop successful postattack sage if aroundsuccfailms==1, absorb(fips i.div_9_all##i.year smale sweapon) cluster(fips)
reghdfe ln_est_pop successful postattack sage if aroundsuccfailms==1, absorb(fips i.div_9_all##i.year smale sweapon) cluster(fips)

******************************************************************************************************************************************
*** TABLE A13
******************************************************************************************************************************************

use "$path\finaldata.dta", clear

*** PANEL A
g interaction=.
replace interaction=successful*schoolshooting
reghdfe ln_emp_pop successful schoolshooting interaction if aroundms==1, absorb(fips i.div_9_all##i.year) cluster(fips)
replace interaction=successful*workshooting
reghdfe ln_emp_pop successful workshooting interaction if aroundms==1, absorb(fips i.div_9_all##i.year) cluster(fips)
replace interaction=successful*publicshooting
reghdfe ln_emp_pop successful publicshooting interaction if aroundms==1, absorb(fips i.div_9_all##i.year) cluster(fips)
replace interaction=successful*nonpublicshooting
reghdfe ln_emp_pop successful nonpublicshooting interaction if aroundms==1, absorb(fips i.div_9_all##i.year) cluster(fips)

*** PANEL B
replace interaction=successful*schoolshooting
reghdfe ln_real_qp1_pop successful schoolshooting interaction if aroundms==1, absorb(fips i.div_9_all##i.year) cluster(fips)
replace interaction=successful*workshooting
reghdfe ln_real_qp1_pop successful workshooting interaction if aroundms==1, absorb(fips i.div_9_all##i.year) cluster(fips)
replace interaction=successful*publicshooting
reghdfe ln_real_qp1_pop successful publicshooting interaction if aroundms==1, absorb(fips i.div_9_all##i.year) cluster(fips)
replace interaction=successful*nonpublicshooting
reghdfe ln_real_qp1_pop successful nonpublicshooting interaction if aroundms==1, absorb(fips i.div_9_all##i.year) cluster(fips)

*** PANEL C
replace interaction=successful*schoolshooting
reghdfe ln_est_pop successful schoolshooting interaction if aroundms==1, absorb(fips i.div_9_all##i.year) cluster(fips)
replace interaction=successful*workshooting
reghdfe ln_est_pop successful workshooting interaction if aroundms==1, absorb(fips i.div_9_all##i.year) cluster(fips)
replace interaction=successful*publicshooting
reghdfe ln_est_pop successful publicshooting interaction if aroundms==1, absorb(fips i.div_9_all##i.year) cluster(fips)
replace interaction=successful*nonpublicshooting
reghdfe ln_est_pop successful nonpublicshooting interaction if aroundms==1, absorb(fips i.div_9_all##i.year) cluster(fips)

*********************************************************************************
*** TABLE A14
*********************************************************************************

use "$path\finaldata.dta", clear

*** PANEL A
g pop15t24 = pop15t19+pop20t24
g pop15t65 = pop*(pop15t24+pop25t34+pop35t44+pop45t54+pop55t64)
g lnpop15t65 = 100*log(1+pop15t65)
g lnpop=100*ln(pop)
replace pop65p=pop*pop65p
g lnpop65p = 100*log(1+pop65p)
g popl15 = pop-pop65p-pop15t65
g lnpopl15 = 100*log(1+popl15)
drop if ln_emp_pop==. & ln_real_qp1_pop==. & ln_est_pop==.

reghdfe lnpop successful if aroundms==1, absorb(fips year) cluster(fips)
reghdfe lnpop15t65 successful if aroundms==1, absorb(fips year) cluster(fips)
reghdfe lnpopl15 successful if aroundms==1, absorb(fips year) cluster(fips)
reghdfe lnpop65p successful if aroundms==1, absorb(fips year) cluster(fips)

*** PANEL B
use "$path\finaldata.dta", clear

reghdfe samehouse successful, absorb(fips year) cluster(fips)
reghdfe moved_outside_county successful, absorb(fips year) cluster(fips)
reghdfe moved_outside_state successful, absorb(fips year) cluster(fips)
reghdfe moved_outside_us successful, absorb(fips year) cluster(fips)

*********************************************************************************
*** TABLE A15
*********************************************************************************

use "$path\media_data2.dta", clear

xi: ivreg2 chln_emp_pop (media = nd2) Total_city_minus_ms i.state_numeric i.divisionyear succdeaths, cluster(fips) first 
xi: ivreg2 chln_real_qp1_pop (media_min = nd2) Total_city_minus_ms i.state_numeric i.divisionyear succdeaths, cluster(fips) first 



******************************************************************************************************************************************
*** FIGURE A1
******************************************************************************************************************************************
use "$path\finaldata.dta", clear

g beta_succ_emp = .
g se_succ_emp = .

reghdfe ln_emp_pop successful14-successful17 successful19-successful26 postattack if aroundsuccfailms==1, absorb(fips year) cluster(fips)
foreach var of varlist successful14-successful17 successful19-successful25 {
replace beta_succ_emp = _b[`var'] if `var'==1
replace se_succ_emp = _se[`var'] if `var'==1
}
replace beta_succ_emp = 0 if successful18==1
replace se_succ_emp = 0 if successful18==1
g lbar_succ_emp = beta_succ_emp - 1.64*se_succ_emp
g ubar_succ_emp = beta_succ_emp + 1.64*se_succ_emp

g beta_succ_ern = .
g se_succ_ern = .

reghdfe ln_real_qp1_pop successful14-successful17 successful19-successful26 postattack if aroundsuccfailms==1, absorb(fips year) cluster(fips)
foreach var of varlist successful14-successful17 successful19-successful25 {
replace beta_succ_ern = _b[`var'] if `var'==1
replace se_succ_ern = _se[`var'] if `var'==1
}
replace beta_succ_ern = 0 if successful18==1
replace se_succ_ern = 0 if successful18==1
g lbar_succ_ern = beta_succ_ern - 1.64*se_succ_ern
g ubar_succ_ern = beta_succ_ern + 1.64*se_succ_ern

g beta_succ_est = .
g se_succ_est = .

reghdfe ln_est_pop successful14-successful17 successful19-successful26 postattack if aroundsuccfailms==1, absorb(fips year) cluster(fips)
foreach var of varlist successful14-successful17 successful19-successful25 {
replace beta_succ_est = _b[`var'] if `var'==1
replace se_succ_est = _se[`var'] if `var'==1
}
replace beta_succ_est = 0 if successful18==1
replace se_succ_est = 0 if successful18==1
g lbar_succ_est = beta_succ_est - 1.64*se_succ_est
g ubar_succ_est = beta_succ_est + 1.64*se_succ_est

keep if yearsince>=-5 & yearsince<=5
duplicates drop yearsince, force
foreach var of varlist lbar* ubar* {
replace `var'=0 if `var'==.
}

twoway (connected beta_succ_emp yearsince) (line lbar_succ_emp yearsince, lpattern(dash) lcolor(black)) ///
(line ubar_succ_emp yearsince, lpattern(dash) lcolor(black) yline(0, lcolor(black))), ///
xtitle("Years since Mass Shooting") ytitle("Impact of Successful MS on Employment") ///
ylabel(-5(1)2) xlabel(-5(1)5) legend(off)

twoway (connected beta_succ_ern yearsince) (line lbar_succ_ern yearsince, lpattern(dash) lcolor(black)) ///
(line ubar_succ_ern yearsince, lpattern(dash) lcolor(black) yline(0, lcolor(black))), ///
xtitle("Years since Mass Shooting") ytitle("Impact of Successful MS on Earnings") ///
ylabel(-5(1)2) xlabel(-5(1)5) legend(off)

twoway (connected beta_succ_est yearsince) (line lbar_succ_est yearsince, lpattern(dash) lcolor(black)) ///
(line ubar_succ_est yearsince, lpattern(dash) lcolor(black) yline(0, lcolor(black))), ///
xtitle("Years since Mass Shooting") ytitle("Impact of Successful MS on Establishments") ///
ylabel(-5(1)2) xlabel(-5(1)4) legend(off)

******************************************************************************************************************************************
*** FIGURE A2
******************************************************************************************************************************************
use "$path\finaldata.dta", clear
sort fips year

g beta_succ_emp = .
g se_succ_emp = .

reghdfe ln_emp_pop successful13-successful18 successful20-successful25 if aroundmsneighbor==1, absorb(fips i.div_9_all##i.year) cluster(fips)
foreach var of varlist successful14-successful18 successful20-successful25 {
replace beta_succ_emp = _b[`var'] if `var'==1
replace se_succ_emp = _se[`var'] if `var'==1
}
replace beta_succ_emp = 0 if successful19==1
replace se_succ_emp = 0 if successful19==1
g lbar_succ_emp = beta_succ_emp - 1.64*se_succ_emp
g ubar_succ_emp = beta_succ_emp + 1.64*se_succ_emp

g beta_succ_ern = .
g se_succ_ern = .

reghdfe ln_real_qp1_pop successful13-successful18 successful20-successful25 if aroundmsneighbor==1, absorb(fips i.div_9_all##i.year) cluster(fips)
foreach var of varlist successful14-successful18 successful20-successful25 {
replace beta_succ_ern = _b[`var'] if `var'==1
replace se_succ_ern = _se[`var'] if `var'==1
}
replace beta_succ_ern = 0 if successful19==1
replace se_succ_ern = 0 if successful19==1
g lbar_succ_ern = beta_succ_ern - 1.64*se_succ_ern
g ubar_succ_ern = beta_succ_ern + 1.64*se_succ_ern

g beta_succ_est = .
g se_succ_est = .

reghdfe ln_est_pop successful13-successful18 successful20-successful25 if aroundmsneighbor==1, absorb(fips i.div_9_all##i.year) cluster(fips)
foreach var of varlist successful14-successful18 successful20-successful25 {
replace beta_succ_est = _b[`var'] if `var'==1
replace se_succ_est = _se[`var'] if `var'==1
}
replace beta_succ_est = 0 if successful19==1
replace se_succ_est = 0 if successful19==1
g lbar_succ_est = beta_succ_est - 1.64*se_succ_est
g ubar_succ_est = beta_succ_est + 1.64*se_succ_est

keep if yearsince>=-5 & yearsince<=5
duplicates drop yearsince, force
foreach var of varlist lbar* ubar* {
replace `var'=0 if `var'==.
}

twoway (connected beta_succ_emp yearsince) (line lbar_succ_emp yearsince, lpattern(dash) lcolor(black)) ///
(line ubar_succ_emp yearsince, lpattern(dash) lcolor(black) yline(0, lcolor(black))), ///
xtitle("Years since Mass Shooting") ytitle("Impact of Successful MS on Employment") ///
ylabel(-5(1)2) xlabel(-5(1)5) legend(off)

twoway (connected beta_succ_ern yearsince) (line lbar_succ_ern yearsince, lpattern(dash) lcolor(black)) ///
(line ubar_succ_ern yearsince, lpattern(dash) lcolor(black) yline(0, lcolor(black))), ///
xtitle("Years since Mass Shooting") ytitle("Impact of Successful MS on Earnings") ///
ylabel(-5(1)2) xlabel(-5(1)5) legend(off)

twoway (connected beta_succ_est yearsince) (line lbar_succ_est yearsince, lpattern(dash) lcolor(black)) ///
(line ubar_succ_est yearsince, lpattern(dash) lcolor(black) yline(0, lcolor(black))), ///
xtitle("Years since Mass Shooting") ytitle("Impact of Successful MS on Establishments") ///
ylabel(-5(1)2) xlabel(-5(1)4) legend(off)

******************************************************************************************************************************************
*** FIGURE A3
******************************************************************************************************************************************
use "$path\finaldata.dta", clear

sort fips year
foreach var of varlist ln_emp_pop ln_real_qp1_pop ln_est_pop ln_ca1_pop_1 lnVIOLENT lnPROPERTY male hs racediv poverty gini lnhom_gun lnsui_gun {
g lag`var'=L1.`var'
}

global controls = "lagln_emp_pop lagln_real_qp1_pop lagln_est_pop laglnVIOLENT laglnPROPERTY lagmale laghs lagracediv lagpoverty laggini laglnhom_gun laglnsui_gun capital_state coastal_county large_hub medium_hub"

***
probit ms lagln_ca1_pop_1 $controls i.year, cluster(fips)
predict p_score, pr
sort p_score year lagln_emp_pop lagln_real_qp1_pop lagln_est_pop fips

g matchms = 0
replace matchms = 1 if ms[_n+1]==1 & mse!=1 
replace matchms = 1 if ms[_n+2]==1 & mse!=1 
replace matchms = 1 if ms[_n+3]==1 & mse!=1 
replace matchms = 1 if ms[_n+4]==1 & mse!=1 
replace matchms = 1 if ms[_n+5]==1 & mse!=1 

sort fips year
by fips: egen matchmse = max(matchms)

g beta_succ_emp = .
g se_succ_emp = .

reghdfe ln_emp_pop successful13-successful18 successful20-successful25 if matchmse==1 | mse==1, absorb(fips year) cluster(fips)
foreach var of varlist successful14-successful18 successful20-successful25 {
replace beta_succ_emp = _b[`var'] if `var'==1
replace se_succ_emp = _se[`var'] if `var'==1
}
replace beta_succ_emp = 0 if successful19==1
replace se_succ_emp = 0 if successful19==1
g lbar_succ_emp = beta_succ_emp - 1.64*se_succ_emp
g ubar_succ_emp = beta_succ_emp + 1.64*se_succ_emp

g beta_succ_ern = .
g se_succ_ern = .

reghdfe ln_real_qp1_pop successful13-successful18 successful20-successful25 if matchmse==1 | mse==1, absorb(fips year) cluster(fips)
foreach var of varlist successful14-successful18 successful20-successful25 {
replace beta_succ_ern = _b[`var'] if `var'==1
replace se_succ_ern = _se[`var'] if `var'==1
}
replace beta_succ_ern = 0 if successful19==1
replace se_succ_ern = 0 if successful19==1
g lbar_succ_ern = beta_succ_ern - 1.64*se_succ_ern
g ubar_succ_ern = beta_succ_ern + 1.64*se_succ_ern

g beta_succ_est = .
g se_succ_est = .

reghdfe ln_est_pop successful13-successful18 successful20-successful25 if matchmse==1 | mse==1, absorb(fips year) cluster(fips)
foreach var of varlist successful14-successful18 successful20-successful25 {
replace beta_succ_est = _b[`var'] if `var'==1
replace se_succ_est = _se[`var'] if `var'==1
}
replace beta_succ_est = 0 if successful19==1
replace se_succ_est = 0 if successful19==1
g lbar_succ_est = beta_succ_est - 1.64*se_succ_est
g ubar_succ_est = beta_succ_est + 1.64*se_succ_est

keep if yearsince>=-5 & yearsince<=5
duplicates drop yearsince, force
foreach var of varlist lbar* ubar* {
replace `var'=0 if `var'==.
}

twoway (connected beta_succ_emp yearsince) (line lbar_succ_emp yearsince, lpattern(dash) lcolor(black)) ///
(line ubar_succ_emp yearsince, lpattern(dash) lcolor(black) yline(0, lcolor(black))), ///
xtitle("Years since Mass Shooting") ytitle("Impact of Successful MS on Employment") ///
ylabel(-5(1)2) xlabel(-5(1)5) legend(off)

twoway (connected beta_succ_ern yearsince) (line lbar_succ_ern yearsince, lpattern(dash) lcolor(black)) ///
(line ubar_succ_ern yearsince, lpattern(dash) lcolor(black) yline(0, lcolor(black))), ///
xtitle("Years since Mass Shooting") ytitle("Impact of Successful MS on Earnings") ///
ylabel(-5(1)2) xlabel(-5(1)5) legend(off)

twoway (connected beta_succ_est yearsince) (line lbar_succ_est yearsince, lpattern(dash) lcolor(black)) ///
(line ubar_succ_est yearsince, lpattern(dash) lcolor(black) yline(0, lcolor(black))), ///
xtitle("Years since Mass Shooting") ytitle("Impact of Successful MS on Establishments") ///
ylabel(-5(1)2) xlabel(-5(1)4) legend(off)


******************************************************************************************************************************************
*** FIGURE A4
******************************************************************************************************************************************

use "$path\finaldata.dta", clear
merge 1:1 fips year using "$path\ms2.dta" 

sort fips year
by fips: egen ms2n = sum(ms2)
by fips: egen ms2e = max(ms2)

by fips: egen minyr2 = min(year) if ms2==1
by fips: egen minyr = min(minyr2) if ms2e==1
by fips: g yrms2 = year if ms2==1
by fips: replace yrms2 = yrms2[_n-1] if yrms==.
by fips: replace yrms2 = minyr if yrms2==. & ms2n!=0 & year<=yrms2
g yearsince2 = year - yrms2
g post2 = (yearsince2>=0 & yearsince2!=.)
drop minyr2 minyr
g mspost2 = ms2e*post2

sort fips year
g aroundms2 = 0
forvalues j=-6(1)6{
by fips: replace aroundms2=1 if ms2[_n+`j']==1
}
replace aroundms2=1 if ms2e==0

g definition=2
foreach x in "b" "lb" "ub" {
g emp_`x' = .
g earn_`x' = .
g est_`x' = .
}

reghdfe ln_emp_pop mspost2 if aroundms2==1, absorb(fips i.div_9_all##i.year) cluster(fips)
lincom mspost2
replace emp_b = `r(estimate)'
replace emp_lb = `r(lb)'
replace emp_ub = `r(ub)'

reghdfe ln_real_qp1_pop mspost2 if aroundms2==1, absorb(fips i.div_9_all##i.year) cluster(fips state_all)
lincom mspost2
replace earn_b = `r(estimate)'
replace earn_lb = `r(lb)'
replace earn_ub = `r(ub)'

reghdfe ln_est_pop mspost2 if aroundms2==1, absorb(fips i.div_9_all##i.year) cluster(fips)
lincom mspost2
replace est_b = `r(estimate)'
replace est_lb = `r(lb)'
replace est_ub = `r(ub)'

keep definition *_b *_ub *_lb
compress
save "$path\results_death_cutoff2.dta", replace 

***
use "$path\finaldata.dta", clear
merge 1:1 fips year using "$path\ms3.dta" 

sort fips year
by fips: egen ms3n = sum(ms3)
by fips: egen ms3e = max(ms3)

by fips: egen minyr2 = min(year) if ms3==1
by fips: egen minyr = min(minyr2) if ms3e==1
by fips: g yrms3 = year if ms3==1
by fips: replace yrms3 = yrms3[_n-1] if yrms3==.
by fips: replace yrms3 = minyr if yrms3==. & ms3n!=0 & year<=yrms3
g yearsince3 = year - yrms3
g post3 = (yearsince3>=0 & yearsince3!=.)
drop minyr2 minyr
g mspost3 = ms3e*post3

sort fips year
g aroundms3 = 0
forvalues j=-6(1)6{
by fips: replace aroundms3=1 if ms3[_n+`j']==1
}
replace aroundms3=1 if ms3e==0

g definition=3
foreach x in "b" "lb" "ub" {
g emp_`x' = .
g earn_`x' = .
g est_`x' = .
}

reghdfe ln_emp_pop mspost3 if aroundms3==1, absorb(fips i.div_9_all##i.year) cluster(fips)
lincom mspost3
replace emp_b = `r(estimate)'
replace emp_lb = `r(lb)'
replace emp_ub = `r(ub)'

reghdfe ln_real_qp1_pop mspost3 if aroundms3==1, absorb(fips i.div_9_all##i.year) cluster(fips state_all)
lincom mspost3
replace earn_b = `r(estimate)'
replace earn_lb = `r(lb)'
replace earn_ub = `r(ub)'

reghdfe ln_est_pop mspost3 if aroundms3==1, absorb(fips i.div_9_all##i.year) cluster(fips)
lincom mspost3
replace est_b = `r(estimate)'
replace est_lb = `r(lb)'
replace est_ub = `r(ub)'

keep definition *_b *_ub *_lb
compress
save "$path\results_death_cutoff3.dta", replace 


*** 
use "$path\finaldata.dta", clear

drop if succdeaths==4
sort fips year
g aroundms2 = 0
forvalues j=-6(1)6{
by fips: replace aroundms2=1 if ms[_n+`j']==1
}
replace aroundms2=1 if mse==0


g definition=5
foreach x in "b" "lb" "ub" {
g emp_`x' = .
g earn_`x' = .
g est_`x' = .
}

reghdfe ln_emp_pop mspost if aroundms2==1, absorb(fips i.div_9_all##i.year) cluster(fips)
lincom mspost
replace emp_b = `r(estimate)'
replace emp_lb = `r(lb)'
replace emp_ub = `r(ub)'

reghdfe ln_real_qp1_pop mspost if aroundms2==1, absorb(fips i.div_9_all##i.year) cluster(fips state_all)
lincom mspost
replace earn_b = `r(estimate)'
replace earn_lb = `r(lb)'
replace earn_ub = `r(ub)'

reghdfe ln_est_pop mspost if aroundms2==1, absorb(fips i.div_9_all##i.year) cluster(fips)
lincom mspost
replace est_b = `r(estimate)'
replace est_lb = `r(lb)'
replace est_ub = `r(ub)'

keep definition *_b *_ub *_lb
compress
save "$path\results_death_cutoff5.dta", replace 

*** 
use "$path\finaldata.dta", clear

drop if succdeaths==4 | succdeaths==5
sort fips year
g aroundms2 = 0
forvalues j=-6(1)6{
by fips: replace aroundms2=1 if ms[_n+`j']==1
}
replace aroundms2=1 if mse==0


g definition=6
foreach x in "b" "lb" "ub" {
g emp_`x' = .
g earn_`x' = .
g est_`x' = .
}

reghdfe ln_emp_pop mspost if aroundms2==1, absorb(fips i.div_9_all##i.year) cluster(fips)
lincom mspost
replace emp_b = `r(estimate)'
replace emp_lb = `r(lb)'
replace emp_ub = `r(ub)'

reghdfe ln_real_qp1_pop mspost if aroundms2==1, absorb(fips i.div_9_all##i.year) cluster(fips state_all)
lincom mspost
replace earn_b = `r(estimate)'
replace earn_lb = `r(lb)'
replace earn_ub = `r(ub)'

reghdfe ln_est_pop mspost if aroundms2==1, absorb(fips i.div_9_all##i.year) cluster(fips)
lincom mspost
replace est_b = `r(estimate)'
replace est_lb = `r(lb)'
replace est_ub = `r(ub)'

keep definition *_b *_ub *_lb
compress
save "$path\results_death_cutoff6.dta", replace 

*** PANEL D
use "$path\finaldata.dta", clear

g definition=4
foreach x in "b" "lb" "ub" {
g emp_`x' = .
g earn_`x' = .
g est_`x' = .
}

reghdfe ln_emp_pop mspost if aroundms==1, absorb(fips i.div_9_all##i.year) cluster(fips)
lincom mspost
replace emp_b = `r(estimate)'
replace emp_lb = `r(lb)'
replace emp_ub = `r(ub)'

reghdfe ln_real_qp1_pop mspost if aroundms==1, absorb(fips i.div_9_all##i.year) cluster(fips state_all)
lincom mspost
replace earn_b = `r(estimate)'
replace earn_lb = `r(lb)'
replace earn_ub = `r(ub)'

reghdfe ln_est_pop mspost if aroundms==1, absorb(fips i.div_9_all##i.year) cluster(fips)
lincom mspost
replace est_b = `r(estimate)'
replace est_lb = `r(lb)'
replace est_ub = `r(ub)'

keep definition *_b *_ub *_lb
compress
save "$path\results_death_cutoff4.dta", replace 


use "$path\results_death_cutoff2.dta", replace 
forvalues j=3(1)6{
append using "$path\results_death_cutoff`j'.dta"
}
save "$path\results_death_cutoff.dta", replace 

forvalues j=2(1)6{
capture: erase "$path\results_death_cutoff`j'.dta"
}

use "$path\results_death_cutoff.dta", clear

twoway (scatter emp_b definition) (rcap emp_lb emp_ub definition, yline(0, lpattern(dash) lcolor(black))),  ///
xtitle("Definition of Mass Shooting (Number of Deaths)") ytitle("Impact of MS on Employment") ///
ylabel(-4(1)2) xlabel(2(1)6) legend(off) 

twoway (scatter earn_b definition) (rcap earn_lb earn_ub definition, yline(0, lpattern(dash) lcolor(black))),  ///
xtitle("Definition of Mass Shooting (Number of Deaths)") ytitle("Impact of MS on Earnings") ///
ylabel(-4(1)2) xlabel(2(1)6) legend(off) 

twoway (scatter est_b definition) (rcap est_lb est_ub definition, yline(0, lpattern(dash) lcolor(black))),  ///
xtitle("Definition of Mass Shooting (Number of Deaths)") ytitle("Impact of MS on Establishments") ///
ylabel(-4(1)2) xlabel(2(1)6) legend(off) 

******************************************************************************************************************************************
*** FIGURE A5
******************************************************************************************************************************************

use "$path\FinalData\BRFSS.dta", clear

g timesinceyear = floor(timesince/12)

tab timesinceyear, gen(successfuly)

g samplees = 0
forvalues j=8(1)17 {
replace successfuly`j'=0 if successfuly`j'==.
replace samplees = 1 if successfuly`j'==1
}
replace samplees = 1 if mse==0

g beta_ment = .
g se_ment = .

reghdfe menthlth successfuly8-successfuly12 successfuly14-successfuly17 elem_school atthigh_school high_school att_coll married divorced widowed separated couple male age age_2 [pw=_finalwt] ///
if samplees==1 & retired!=1 & student!=1 & disabled!=1 & homemaker!=1 & age<65, cluster(fips) absorb(fips year)

foreach var of varlist successfuly8-successfuly12 successfuly14-successfuly17 {
replace beta_ment = _b[`var'] if `var'==1
replace se_ment = _se[`var'] if `var'==1
}
replace beta_ment = 0 if successfuly13==1
replace se_ment = 0 if successfuly13==1
g lbar_ment = beta_ment - 1.64*se_ment
g ubar_ment = beta_ment + 1.64*se_ment

g beta_psy = .
g se_psy = .

reghdfe physhlth successfuly8-successfuly12 successfuly14-successfuly17 elem_school atthigh_school high_school att_coll married divorced widowed separated couple male age age_2 [pw=_finalwt] ///
if samplees==1 & retired!=1 & student!=1 & disabled!=1 & homemaker!=1 & age<65, cluster(fips) absorb(fips year)

foreach var of varlist successfuly8-successfuly12 successfuly14-successfuly17 {
replace beta_psy = _b[`var'] if `var'==1
replace se_psy = _se[`var'] if `var'==1
}
replace beta_psy = 0 if successfuly13==1
replace se_psy = 0 if successfuly13==1
g lbar_psy = beta_psy - 1.64*se_psy
g ubar_psy = beta_psy + 1.64*se_psy

g beta_ovrl = .
g se_ovrl = .

reghdfe poorhlth successfuly8-successfuly12 successfuly14-successfuly17 elem_school atthigh_school high_school att_coll married divorced widowed separated couple male age age_2 [pw=_finalwt] ///
if samplees==1 & retired!=1 & student!=1 & disabled!=1 & homemaker!=1 & age<65, cluster(fips) absorb(fips year)

foreach var of varlist successfuly8-successfuly12 successfuly14-successfuly17 {
replace beta_ovrl = _b[`var'] if `var'==1
replace se_ovrl = _se[`var'] if `var'==1
}
replace beta_ovrl = 0 if successfuly13==1
replace se_ovrl = 0 if successfuly13==1
g lbar_ovrl = beta_ovrl - 1.64*se_ovrl
g ubar_ovrl = beta_ovrl + 1.64*se_ovrl

duplicates drop timesinceyear, force
keep if timesinceyear>=-5 & timesinceyear<=4

twoway (connected beta_ment timesinceyear) (line lbar_ment timesinceyear, lpattern(dash) lcolor(black)) ///
(line ubar_ment timesinceyear, lpattern(dash) lcolor(black) yline(0, lcolor(black))), ///
xtitle("Years since Mass Shooting") ytitle("Impact of Successful MS on Mental Health") ///
ylabel(-1.5(0.5)1.5) xlabel(-5(1)4) legend(off)

twoway (connected beta_psy timesinceyear) (line lbar_psy timesinceyear, lpattern(dash) lcolor(black)) ///
(line ubar_psy timesinceyear, lpattern(dash) lcolor(black) yline(0, lcolor(black))), ///
xtitle("Years since Mass Shooting") ytitle("Impact of Successful MS on Physical  Health") ///
ylabel(-1.5(0.5)1.5) xlabel(-5(1)4) legend(off)

twoway (connected beta_ovrl timesinceyear) (line lbar_ovrl timesinceyear, lpattern(dash) lcolor(black)) ///
(line ubar_ovrl timesinceyear, lpattern(dash) lcolor(black) yline(0, lcolor(black))), ///
xtitle("Years since Mass Shooting") ytitle("Impact of Successful MS on Overall Health") ///
ylabel(-1.5(0.5)1.5) xlabel(-5(1)4) legend(off)





